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Abstract. We adapt the CRT approach for computing Hilbert class 
polynomials to handle a wide range of class invariants. For suitable dis- 
criminants D, this improves its performance by a large constant factor, 
more than 200 in the most favourable circumstances. This has enabled 
record-breaking constructions of elliptic curves via the CM method, in- 
cluding examples with \D\ > 10^^. 

1 Introduction 

Every ordinary elliptic curve E over a finite field has complex multiplication 
by an imaginary quadratic order O, by which we mean that the endomorphism 
ring End(£') is isomorphic to O. The Deuring lifting theorem implies that E 
is the reduction of an elliptic curve E/C that also has complex multiplication 
by O. Let K denote the fraction field of O. The j-invariant of E is an algebraic 
integer whose minimal polynomial over K is the Hilbert class polynomial H^, 
where D is the discriminant of O. Notably, the polynomial Hd actually lies in 
Z[X], and its splitting field is the ring class field Kq for the order O. 

Conversely, an elliptic curve E /¥q with complex multiplication by O exists 
whenever q satisfies the norm equation Aq — t^ — v^D, with t,?; S Z and i ^ 
modulo the characteristic of Fg. In this case Ho splits completely over Fg, and its 
roots are precisely the j-invariants of the elliptic curves E/Vq that have complex 
multiplication by O. Such a curve has q+l±t points, where t is determined, 
up to a sign, by the norm equation. With a judicious selection of D and q one 
may obtain a curve with prescribed order. This is known as the CM method. 

The main challenge for the CM method is to obtain the polynomial Hd , which 
has degree equal to the class number h{D), and total size 0(|£'|^+'^). There are 
three approaches to computing Hjy, all of which, under reasonable assumptions, 
can achieve a running time of 0{\D\^^'^). These include the complex analytic 
method [12], a p-adic algorithm [9, 7], and an approach based on the Chinese 
Remainder Theorem (CRT) [2]. The first is the most widely used, and it is quite 
efficient; the range of discriminants to which it may be applied is limited not by 
its running time, but by the space required. The polynomial Hjj is already likely 
to exceed available memory when \D\ > 10^, hence one seeks to apply the CM 
method to alternative class polynomials that have smaller coefficients than Ho. 
This makes computations with \D\ > 10^" feasible. 

Recently, a modified version of the CRT approach was proposed that greatly 
reduces the space required for the CM method [30] . Under the Generalised Rie- 
mann Hypothesis (GRH), this algorithm is able to compute Hjj mod P using 



0{\D\^/'^~^'^ logP) space and 0{\D\^~^'^) time. (Here and in the following, all com- 
plexity estimates refer to bit operations.) The reduced space complexity allows 
it to handle much larger discriminants, including examples with \D\ > 10^'^. 

An apparent limitation of the CRT approach is that it depends on some 
specific features of the j-function. As noted in [2] , this potentially precludes it 
from computing class polynomials other than Hjj. The purpose of the present 
article is to show how these obstructions may be overcome, allowing us to apply 
the CRT method to many functions other than j, including two infinite families. 

Subject to suitable constraints on D, we may then compute a class polynomial 
with smaller coefficients than Hu (by a factor of up to 72), and, in certain cases, 
with smaller degree (by a factor of 2). Remarkably, the actual running time with 
the CRT method is typically better than the size difference would suggest. Fewer 
CRT moduli are needed, and we may choose a subset for which the computation 
is substantially faster than on average. 

We start §2 with a brief overview of the CRT method, and then describe a 
new technique to improve its performance, which also turns out to be crucial for 
certain class invariants. After discussing families of invariants in §3, we consider 
CRT-based approaches applicable to the different families and give a general 
algorithm in §4. Computational results and performance data appear in §5. 

2 Hilbert class polynomials via the CRT 

2.1 The Eilgorithm of Belding, Broker, Enge, Lauter and Sutherland 

The basic idea of the CRT-based algorithm for Hilbert class polynomials is to 
compute Hu modulo many small primes and then lift its coefficients by Chi- 
nese remaindering to integers, or to their reductions modulo a large (typically 
prime) integer P, via the explicit CRT [4, Thm. 3.1]. The latter approach suf- 
fices for most applications, and while it does not substantially reduce the running 
time (the same number of small primes is required) , it can be accomplished using 
only 0(|D|i/2+'=iogP) space with the method of [30, §6]. 

For future reference, we summarise the algorithm to compute Hu mod p for 
a prime p that splits completely in the ring class field Kq- Let h = h{D). 

Algorithm 1 (Computing Hd mod p) 

1. Find the j -invariant ji of an elliptic curve E/¥p with End(£^) = O. 

2. Enumerate the other roots j2, - ■ ■, jh of Hd mod p. 

3. Compute Hn{X) mod p = {X - ji) ■ ■ ■ {X - jh). 

The first step is achieved by varying ji (systematically or randomly) over the 
elements of Fp until it corresponds to a suitable curve; details and many practical 
improvements are given in [2, 30]. The third step is a standard building block of 
computer algebra. Our interest lies in Step 2. 
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2.2 Enumerating the roots of Hd mod p 

The key idea in [2] leading to a quasi-hnear complexity is to apply the Galois 
action of C\{0) ~ G&\{Ko/K). The group 01(0) acts on the roots of Hd, 
and when p splits completely in Kq there is a corresponding action on the set 
Ellc)(Fp) = {ji, . . . , j/i} containing the roots of Hd mod p. For an ideal class [a] 
in 01(0) and a j-invariant ji e Ellc)(Fp), let us write [a\ji for the image of ji 
under the Galois action of [a]. Wc then have Elle)(Fp) = {[a]ji : [a] S 01(O)}. 

As in [30, §5], wc use a polycyclic presentation defined by a sequence of 
ideals [i, . . . , with prime norms i\,...,irn whose classes generate 01(O). The 
relative order rk is the least positive integer for which [[^''] e ([Ii], . . . , [Ife-i])- 
We may then uniquely write [a] = [t'l^] • • • [If™], with < < r^. To maximise 
performance, wc use a presentation in which (.i < ••• < with each as 
small as possible subject to > 1. Note that the relative order divides the 
order rik of [ik] in 01(0), but for > 1 we can (and often do) have rk < rik- 

For each G Ellc>(Fp) and each O-ideal [ of prime norm £, the )-invariant 
corresponds to an £-isogenous curve, which we may obtain as a root of ^e{ji,X), 
where S Z[J, J^] is the classical modular polynomial [31, §69]. The polyno- 
mial $^ has the pair of functions [j{z),j{£z)) as roots, and parameterises isoge- 
nics of degree i. 

Fixing an isomorphism End{E) = O, we let tt e O denote the Frobenius 

endomorphism. When the order Z[7r] is maximal at £, the univariate polynomial 
^i{ji,X) G Fp[X] has exactly two roots and when i splits in O, and 
a single root if £ is ramified [25, Prop. 23]. To simplify matters, we assume 

here that Z[7r] is maximal at each £k, but this is not necessary, see [30, §4]. 
Wc may enumerate Ellc)(Fp) = {[a]ji : [a] G ([li], • • • , [U])} via [30, Alg. 1.3]: 

Algorithm 2 (Enumerating Ellcj(Fp) — Step 2 of Algorithm 1) 

1. Let j2 be an arbitrary root of X) in Fp. 

2. Fori from 3 to r^, let ji be the root of ^^^{ji-i, X) / {X — ji_2) in¥p. 

3. If m> 1, then for i from 1 to rm: 

Recursively enumerate the set {[a\ji : [a] G ([Ii], . . . , [(m_i])}. 

In general there are two distinct choices for j2, but either will do. Once j2 is 
chosen, js, . . . , arc determined. The sequence (ji, . . . , jV„) corresponds to a 
path of £m-isogenies; we call this path an £^-thread. 

The choice of j2 in Step 1 may change the order in which Elle)(Fp) is enumer- 
ated. Three of the sixteen possibilities when m — 2, ri — A, and r2 = 3 arc shown 
below; we assume [[|] = [li], and label each vertex [l2]ji by the exponent e. 




Bold edges indicate where a choice was made. Regardless of these choices. 
Algorithm 2 correctly enumerates Ello(Fp) in every case [30, Prop. 5]. 
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2.3 Finding roots with greatest common divisors (gcds) 



The potentially haphazard manner in which Algorithm 2 enumerates Ellc)(Fp) is 
not a problem when computing Hu , but it can complicate matters when we wish 
to compute other class polynomials. We could distinguish the actions of I and T 
using an Elkies kernel polynomial [10], as suggested in [7, §5], however this slows 
down the algorithm significantly. An alternative approach using polynomial gcds 
turns out to be much more efficient, and actually speeds up Algorithm 2, making 
it already a useful improvement when computing Hd. 

We need not distinguish the actions of i and I at this stage, but we wish to 
ensure that our enumeration of Elle)(Fp) makes a consistent choice of direction 
each time it starts an ^-thread. The first ^-thread may be oriented arbitrarily, 
but for each subsequent £-thread (j'i, j2, ■ ■ ■ we apply Lemma 1 below. This 
allows us to "square the corner" by choosing j'2 as the unique common root of 
<^t{X^i'^) and <I>£/(X, j2), where (ji, . . . , jr) is a previously computed ^-thread 
and ji is £'-isogenous to j[. The edge (ji, jj) lies in an ^'-thread that has already 
been computed, for some £' > i. 



Having computed J2, we could compute jg, . . . , as before, but it is usually 
better to continue using gcds, as depicted above. Asymptotically, both root- 
finding and gcd computations are dominated by the 0(^^M(logp)) time it takes 
to instantiate ^)>{X, ji) mod p, but in practice £ is small, and we effectively gain 
a factor of O(logp) by using gcds when i £' . This can substantially reduce the 
running time of Algorithm 2, as may be seen in Table 1 of §5. 

With the gcd approach described above, the total number of root-finding 
operations can be reduced from HfeLi to J2k=i ^k- When m is large, this is a 
big improvement, but it is no help when m = 1, as necessarily occurs when h{D) 
is prime. However, even in this case we can apply gcds by looking for an auxiliary 
ideal I'l, with prime norm £[, for which = [(^]. When n is large, such an l'^ is 
easy to find, and we may choose the best combination of £[ and e available. This 
idea generalises to ^fe-threads, where we seek [['^] e ([li] . . . , [tfe])\([ti] • • • , ['fe-i])- 

Lemma 1. Let 31,32 S Ello(Fp), and let ^1,^2 ^ p be distinct primes with 
A£l£2 < \D\. Then gcd($^;^(ji,X), $^2(^2,^)) has degree at most 1. 

Proof. It follows from [25, Prop. 23] that ji) and ^^^(X, ^2) have at most 

two common roots in the algebraic closure Fp, which in fact lie in Ellc)(Fp). If 
there arc exactly two, then both £1 = (1(1 and £2 = hh split in O, and one of [f (3 
or [JI2 is principal with a non-rational generator. We thus have a norm equation 
4^f^2 = 0? - b'^D with a,b and 6^0, and the lemma follows. 
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3 Class invariants 



Due to the large size of Ho , much effort has been spent seeking smaller generators 
of Kq- For a modular function / and O = Z[t], with r in the upper half plane, 
we call /(t) a class invariant if /(r) e Kq- The class polynomial for / is 

HD[f]iX)= n (^-[a]/(r)). 

Meci(O) 

The contemporary tool for determining class invariants is Shimura's reciprocity 
law; see [28, Th. 4] for a fairly general result. Class invariants arising from many 
different modular functions have been described in the literature; we briefly 
summarise some of the most useful ones. 

Let rj be Dedekind's function, and let = exp{2ni/n). Weber considered 

fM-'?(i) w.^-^/9^(2^) 

^^^"^-W ^'^"^""^^W 

powers of which yield class invariants when (y) ^ — 1, and also 72 = which 
is a class invariant whenever 3 ] D. The Weber functions can be generalised 
[15, 16, 21, 20, 23], and we have the simple and double //-quotients 

n,^(.) = ^^^^^^ = ^{^)^ {^) ^.^^ ^ ^ ^^^^^ 

T]{Z) 



where pi and P2 are primes. Subject to constraints on D, including that no prime 
dividing A'' is inert in O, suitable powers of these functions yield class invariants, 
see [15, 16]. For s = 24/gcd(24, {pi — l){p2 — 1)), the canonical power tu*^ 
is invariant under the Fricke involution W\n ■ z 1-^ for r°(iV), equivalently, 
the Atkin-Lehner involution of level TV, by [17, Thm. 2]. 

The theory of [28] applies to any functions for T^{N), in particular to those of 
prime level N invariant under the Fricke involution, which yield class invariants 
when (^) 7^—1. Atkin developed a method to compute such functions Ajq, which 
are conjectured to have a pole of minimal ordc^r at the unique cusp [10, 26]. These 
are used in the SEA algorithm, and can be found in A4AGMA or Pari/GP. 

The functions above all yield algebraic integers, so Hnlf] G C/f [-'^]- Except 
for tr^ or when gcd(A'^, D) 7^ 1, in which cases additional restrictions may apply, 
one actually has HdU] G '^[X]^ cf. [16, Cor. 3.1]. The (logarithmic) height of 
HdU] = O'iX^ is log max \ai\, which determines the precision needed to com- 
pute the ai. We let Coif) denote the ratio of the heights of HdIJ] and Hnlf]. 

With c(/) = lim|c|^oo coif), we have: 0(72) - 3; c(f) = 72 (when (f ) = 1); 

, 24(7V+1) ^ , ^ 12V'(piP2) ... A^ + 1 

^^•^-^ = 4N^-^ = .(Pi-l)(P2-l) ' '^^^^ = WM' 

where e divides the exponent s defined above, is the order of the pole of An 
at the cusp, and il){pipi) is {p\ + l){pi + 1) when p\ ^ p2, and pi{pi + 1) when 



5 



Pi = p2- Morain observed in [27] that 0(^71) = 36, which is so far the best value 
known when (^) = —1. We conjecture that in fact for all primes TV > 11 with 
iV = 11 mod 60 we have c{An) = 30;^^, and that for N = -I mod 60 we 
have c(j4jv) = 30. This implies that given an arbitrary discriminant D, we can 
always choose N so that An yields class invariants with cd{An) > 30 + o(1). 

When the prime divisors of N are all ramified in K, both IDpj.pj and An 
yield class polynomials that are squares in Z[X], see [11, §1.6] and [18]. Taking 
the square root of such a class polynomial reduces both its degree and its height 
by a factor of 2. For a composite fundamental discriminant D (the most common 
case), this applies to Hd[An] for any prime N \ D.ln the best case, D is divisible 
by 71, and we obtain a class polynomial that is 144 times smaller than Hn- 

3.1 Modulcir polynomials 

Each function f(z) considered above is related to j{z) by a modular polynomial 
VE"/ € Z[i^, J] satisfying ^ f{f{z),j{z)) = 0. For primes £ not dividing the level N, 
we let ^ij denote the minimal polynomial satisfying ^e,f{f{z), f{(-z)) = 0; it is 
a factor of Resj^ (Resj($^(J, Ji),^f{F, J)), ^/(F^, Je)), and as such, an element 
of Z[F,Fe]. Thus generalises the classical modular polynomial ^£ = ^e.j- 

The polynomial ^ij has degree d{£+l) in F and F^, where d divides degj ^f/, 
see [6, §6.8], and 2d divides degj^/ when / is invariant under the Pricke invo- 
lution. In general, d is maximal, and d = 1 is achievable only in the relatively 
few cases where Xo{N), respectively Xq{N), is of genus and, moreover, / is 
a hauptmodul, that is, it generates the function field of the curve. Happily, this 
includes many cases of practical interest. 

The polynomial \1/ / characterises the analytic function / in an algebraic way; 
when d = 1, the polynomials and ^ej algebraically characterise ^-isogenies 
between elliptic curves given by their j-invariants, or by class invariants derived 
from /, respectively. These are key ingredients for the CRT method. 

4 CRT algorithms for class invariants 

To adapt Algorithm 1 to class invariants arising from a modular function /(z) 
other than j{z), we only need to consider Algorithm 2. Our objective is to 
enumerate the roots of Hd [/] mod p for suitable primes p, which we are free 
to choose. This may be done in one of two ways. The most direct approach 
computes an "/-invariant" /i, corresponding to ji, then enumerates f2, ■ ■ ■ , fh 
using the modular polynomials ^e,f- Alternatively, we may enumerate ji,. ■ ■ ,jh 
as before, and from these derive /i, . . . , f^. The latter approach is not as efficient, 
but it applies to a wider range of functions, including two infinite families. 

Several problems arise. First, an elliptic curve E/¥p with CM by O unam- 
biguously defines a j'-invariant ji = j{E), but not the corresponding /i. The /i 
we seek is a root of ^f{X) = modp, but ^/i/ may have other roots, 

which may or may not be class invariants. The same problem occurs for the 
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p-adic lifting algorithm and can be solved generically [6, §6]; we describe some 
more efficient solutions, which arc in part specific to certain types of functions. 

When il>f has multiple roots that arc class invariants, these may be roots 
of distinct class polynomials. We are generally happy to compute any one of 
these, but it is imperative that we compute the reduction of "the same" class 
polynomial [/] modulo each prime p. 

The lemma below helps to address these issues for at least two infinite families 
of functions: the double Tj-quotients ttipi^p^ and the Atkin functions Ajv- 

Lemma 2. Let f he a modular function for T^{N), invariant under the Fricke 

involution W\m , such that f{z) and f {--^) have rational q-expansions. Let the 
imaginary quadratic order O have conductor coprime to N and contain an 
ideal n = {N, ^°+^ ). Let Aq = and tq = and assume that 

gcd{Ao, N) — 1. Then /(tq) is a class invariant, and if /(t) is any of its conju- 
gates under the action of Gal{Ko/ K) we have 

=0 and %(/(T),[n]i(r)) =0. 

Proof. By definition, ^'/(/(z), = 0. Applying the Fricke involution yields 

= ^f{{W\Nf){z),{W\Nj){z)) = % {f{z),j{=f)) = M// The 
value /(to) is a class invariant by [28, Th. 4]. By the same result, we may assume 
that T is the basis quotient of an ideal o = {A, ~^~^^ ) with gcd(^, N) = 1 
and B = Bq mod 2N. Then ^ is the basis quotient of an = {AN, ~-^+^ ). It 
follows that [n]j(T) = j (jf), and replacing z above by r completes the proof. 

If wc arrange the roots of Hd into a graph of n-isogeny cycles corresponding 
to the action of n, the lemma yields a dual graph defined on the roots of Hnlf], 
in which vertices /(r) correspond to edges (jir), [Ti]j(r)). 

In computational terms, /(r) is a root of gcd f{X,j{T)) ^'i' f{X, [n]j{T))). 
Generically, we expect this gcd to have no other roots modulo primes p that split 
completely in Ko- For a finite number of such primes, there may be additional 
roots. We have observed this for p dividing the conductor of the order generated 
by /(r) in the maximal order of Ko- Such primes may either be excluded from 
our CRT computations, or addressed by one of the techniques described in §4.3. 

4.1 Direct enumeration 

When the polynomials j have degree £ + 1 we can apply Algorithm 2 with 
essentially no modification; the only new consideration is that i must not divide 
the level N, but we can exclude such £ when choosing a polycyclic presentation 
for 01(0). When the degree is greater than £+1 the situation is more complex, 
moreover the most efficient algorithms for computing modular polynomials do 
not apply [8, 13], making it difficult to obtain unless £ is very small. Thus 
in practice we do not use j in this case; instead we apply the methods of §4.3 
or §4.4. For the remainder of this subsection and the next we assume that we do 
have polynomials ^ej of degree i + 1 with which to enumerate fi,. . .,fh, and 
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consider how to determine a starting point /i, given the j-invariant ji = j{E) 
of an eUiptic curve E/¥p with CM by O. 

When il>f{X) = ^'/(X, ji) mod p has only one root, our choice of /i is imme- 
diately determined. This is usually not the case, but we may be able to ensure it 
by restricting our choice of p. As an example, for / = 72 with 3 \ D, if we require 
that p = 2 mod 3, then /i is the unique cube root of ji in ¥p. If we addition- 
ally have D = 1 mod 8 and p = 3 mod 4, then the equation 72 = (f^'' — 16)/f* 
uniquely determines the square of the Weber f function, by [8, Lem. 7.3]. To 
treat f itself we need an additional trick described in §4.2. 

The next simplest case occurs when only one of the roots of V'/ is a class 
invariant. This necessarily happens when / is invariant under the Fricke involu- 
tion and all the primes dividing A'' are ramified in O. In the context of Lemma 2, 
each root of Holf] then corresponds to an isolated edge (jXt), [njjXT)) in the 
n-isogcny graph on the roots of Hjj, and we compute /i as the unique root of 
gcd(^'/(X, ji), ^'/(A, [n]ji)) . In this situation n = n, and each /(t) occurs twice 
as a root of Holf]. By using a polycyclic presentation for Cl(C')/([n]) rather 
than C1(C), we enumerate each double root of Hjj[f] mod p just once. 

Even when tpf has multiple roots that are class invariants, it may happen 
that they are all roots of the same class polynomial. This applies to the Atkin 
functions f — A^. When A^ is a split prime, there are two A^-isogenous pairs 
(ji,[n]ji) and {[n]ji,ji) in Elle)(Fp), and under Lemma 2 these correspond to 
roots /i and [fi]/i of V'/- Both are roots of Hoif], and we may choose either. 

The situation is slightly more complicated for the double 77-quotients tVp^ , 
with N = P1P2 composite. If pi = pipi and p2 = P2P2 both split and pi ^ p2, 
then there are four distinct A'-isogenies corresponding to four roots oi ipf- Two 
of these roots are related by the action of [n] = [P1P2]; they belong to the same 
class polynomial, which we choose as Hd [f] mod p. The other two are related 
by [P1P2] and are roots of a different class polynomial. We make an arbitrary 
choice for /i, explicitly compute [n]/i, and then check whether it occurs among 
the other three roots; if not, we correct the initial choice. The techniques of §4.3 
may be used to efficiently determine the action of [n]. 

Listed below arc some of the modular functions / for which the roots of 
Hulf] mod p may be directly enumerated, with sufficient constraints on D and p. 
In each case p splits completely in Kq and D < — 4A'^ has conductor u. 

(1) 72, with 3 t -D and p = 2 mod 3; 

(2) f, with D=l mod 8,3\D, and p = 11 mod 12; 

(3) ro%, for N e {3, 5, 7, 13} and s = 24/ gcd(24, N - 1), with N\D and N\u; 

(4) rof , with 3] D,5\ D. and 5 f u; 

(5) An, for A^ e {3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 41, 47, 59, 71}, with (f) ^ -1 
and N ]u. 

(6) ro;^^p^, for {puP2) e {(2, 3), (2, 5), (2, 7), (2, 13), (3, 5), (3, 7), (3, 13), (5, 7)} 
and s = 24/gcd(24, (pi - l)(p2 - 1)), with ^ -1 and Pi,P2 

(7) tt)|_3 with (f ) = 1 and 3 1 M. 
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4.2 The trace trick 



In §4.1 we were able to treat the square of the Weber f function but not f itself. 
To remedy this, we generalise a method suggested to us by Reinier Broker. 

Wc consider the situation where there are two modular functions / and /' 
that are roots of j(2;)), both of which yield class invariants for O, and 

we wish to apply the direct enumeration approach. We assume that p is chosen 
so that tpf{X) = 'i>f{X, ji) mod p has exactly two roots, and depending on 
which root we take as /i, we may compute the reduction of either 
or (X) modulo p. In the case of Weber f, we have /' = — /, and Holf] 

differs from Hu[f] only in the sign of every other coefficient. 

Consider a fixed coefficient of (X) = ^a^X*; most of the time, 

the trace t = —ah-i = fi + ■ ■ ■ + fh will do (if /' = — /, we need to use 
with i ^ h mod 2). The two roots /i and /{ lead to two possibilities t and 
modulo p. However, the elementary symmetric functions Ti = t + t' and T2 = tt' 
are unambiguous modulo p. Computing these modulo many primes p yields Ti 
and T2 as integers (via the CRT), from which t and t' arc obtained as roots of 
the quadratic equation X^ — TiX + T2. If these arc different, we arbitrarily pick 
one of them, which, going back, determines the set of conjugates {fi, ■ ■ ■ , fh} or 
{/(, . . . , f'f^} to take modulo each of the primes p \ t — t' . In the unlikely event 
that they are the same (the suspicion t = t' being confirmed after, say, looking 
at the second prime), we need to switch to a different coefficient Oj. 

If / and /' differ by a simple transformation (such as /' = — /), the second 
set of conjugates and the value t' are obtained essentially for free. As a special 
case, when h is odd and the class invariants are units (as with Weber f), we can 
simply fix i = oq = 1, and need not compute Ti = and T2 = —1. 

The key point is that the number of primes p we use to determine t is much 
less than the number of primes we use to compute Hnlf]- Asymptotically, the 
logarithmic height of the trace is smaller than the height bound we use for 
Hoif] by a factor quasi-linear in log|£)|, under the GRH. In practical terms, 
determining t typically requires less than one tenth of the primes used to compute 
Hjylf], and these computations can be combined. 

The approach described above generalises immediately to more than two 
roots, but this case does not occur for the functions we examine. Unfortunately 
it can be used only in conjunction with the direct enumeration approach of §4.1; 
otherwise we would have to consistently distinguish not only between /i and /( , 
but also between fi and f- for i = 2, . . . , /i. 

4.3 Enumeration via the Fricke involution 

For functions / to which Lemma 2 applies, we can readily obtain the roots 
of -&£)[/] modj3 without using the polynomials ^ej- We instead enumerate the 
roots of Hjj mod p (using the polynomials ^e), and arrange them into a graph G 
of n-isogeny cycles, where n is the ideal of norm appearing in Lemma 2. We 
then obtain roots of Hoif] mod p by computing gcd(^'/(X, jj), ^'/(X, [n]ji)) for 
each edge {ji, [n]ji) in G. 
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The graph G is composed of h/n cycles of length n, where n is the order of 
[n] in C1(C). We assume that the O-ideals of norm A'' are all non-principal and 
inequivalent (by requiring \D\> AN^ if needed). When every prime dividing A'' 
is ramified in O we have n = 2; as noted in §4.1, every root of Hjjif] then occurs 
with multiplicity 2, and we may compute the square-root of Holf] by taking 
each root just once. Otherwise we have n > 2. 

Let [li], . . . , [Ijn] be a polycyclic presentation for 01(0) with relative orders 
ri,...,rm, as in §2.2. For k from 1 to m let us fix 1^ = [£k, ~^^+^ ) with 
Bk > 0. To each vector e = (ei, . . . , em) with < < r^, we associate a unique 
root je enumerated by Algorithm 2, corresponding to the path taken from ji 
to je, where Cfc counts steps taken along an ^fc-thread. For o = (0, . . . , 0) we have 
jo = ji , and in general 

with cTfc = ±1. Using the method of §2.3 to consistently orient the £fc-threads 
ensures that each ak depends only on the orientation of the first £fe-thread. 

To compute the graph G we must determine the signs CTfe. For those [Ik] of 
order 2, we let Ufe = 1. We additionally fix cta; = 1 for the least k = ko (if any) for 
which [Ik] has order greater than 2, since we need not distinguish the actions of n 
and fi. It suflaces to show how to determine ak, given that we know cti, . . . , (Tk-i- 
We may assume [Ikg] and [Ik] both have order greater than 2, with ko < k < m. 

Let I be an auxiliary ideal of prime norm £ such that [[] = [ab] = [11^ ' ' ' 'fc*]) 
with < ei < ri, where b = l^*, and [a] and [b] have order greater than 2. Our 
assumptions guarantee that such an [ exists, by the Cebotarev density theorem, 
and under the GRH, £ is relatively small [1] . The fact that [a] and [b] have order 
greater than 2 ensures that [ab] is distinct from [I] and its inverse. It follows that 
Ufe = 1 if and only if ^e{jo, je) = 0, where e = (ei, . . . , Cfc, 0, . . . , 0). 

Having determined the ak, we compute the unique vector v = {vi, . . . , w,„) for 
which [n] = [[^^"^ • • • l^""']. We then have [n]jo = jv, yielding the edge {jo,jv) 
of G. In general, we obtain the vector corresponding to [n]je by computing e + v 
and using relations [[^I*] — [l^^ ■ ■ ■ ik'^T^] to reduce the result, cf. [30, §5]. 

This method may be used with any function / satisfying Lemma 2, and in 
particular it applies to two infinite families of functions: 

(8) Aiv, for TV > 2 prime, with (^) ^ -1 and N \ u. 

(9) "'pi.ps' for Pi'P2 primes not both 2, with (^), (^) ^ -1 and pi,p2 j"- 

As above, u denotes the conductor of D < —AN"^. 

As noted earlier, for certain primes p we may have difficulty computing the 
edges of G when gcd(\l//(X, jj), \1//(X, [n]jj)) has more than one root in F^. 
While we need not use such primes, it is often easy to determine the correct 
root. Here we give two heuristic techniques for doing so. 

The first applies when is prime, as with the Atkin functions. In this case 
problems can arise when Hd [/] has repeated roots modulo p. By Rummer's cri- 
terion, this can happen only when p divides the discriminant of Hd [/] , and even 
then, a repeated root xx is only actually a problem when it corresponds to two al- 
ternating edges in G, say (ji, J2) and {j3,j4), with the edge (j2, Js) between them. 
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In this scenario we will get two roots xi and X2 of gcd(\E'/(X, j2), ^'/(^, Js))- 
But if we already know that Xi corresponds to (ji, ^2), we can unambiguously 
choose X2- In each of the A''-isogeny cycles of G, it is enough to find a single edge 
that yields a unique root. If no such edge exists, then every edge must yield the 
same two roots Xi and X2, and we count each with multiplicity n/2. 

The second technique applies when the roots of Hjj [f] are units, as with the 
double //-quotients [16, Thm. 3.3]. The product of the roots is then ±1. Assuming 
that the number of edges in G for which multiple roots arise is small (it is usually 
zero, and rarely more than one or two), we simply test all the possible choices 
of roots and see which yield ±1. If only one combination works, then the correct 
choices are determined. This is not guaranteed to happen, but in practice it 
almost always does. 

4.4 A general algorithm 

We now briefly consider the case of an arbitrary modular function / of level N, 
and sketch a general algorithm to compute -&£)[/] with the CRT method. 

Let us assume that /(r) is a class invariant, and let D be the discriminant 
and u the conductor of the order O = [l,r]. The roots of ^'/(X, j(T)) e ife)[X] 
lie in the ray class field of conductor uN over K, and some number n of these, 
including /(r), actually lie in the ring class field Kq- We may determine n using 
the method described in [6, §6.4], which computes the action of {O/NO)* /O* on 
the roots of "ii f(X,j{T)). We note that the complexity of this task is essentially 
fixed as a function of ]£>]. 

Having determined n, we use Algorithm 2 to enumerate the roots ji , . . . , jh of 
Hd mod p as usual, but if for any ji we find that ^f{X, ji) mod p does not have 
exactly n roots f^^^ , ■ ■ ■ , fl^^ , we exclude the prime p from our computations. The 
number of such p is finite and may be bounded in terms of the discriminants 
of the polynomials ^ f{X,a) as a ranges over the roots of Holf]- We then 

compute the polynomial H{X) = HiLi 0^=1 ~ fi^^) °f degree nh in Fp[X]. 
After doing this for sufficiently many primes p, we can lift the coefficients by 
Chinese remaindering to the integers. The resulting fl" is a product of n distinct 
class polynomials, all of which may be obtained by factoring H in Z[X]. Under 
suitable heuristic assumptions (including the GRH), the total time to compute 
Ho[f] is quasi-linear in \D\, including the time to factor H. 

This approach is practically efficient only when n is small, but then it can 
be quite useful. A notable example is the modular function g for which 

*g(X, J) = (A:^^ _ _ 27)3 _ jx^». 

This function was originally proposed by Atkin, and is closely related to certain 
class invariants of Ramanujan [3, Thm. 4.1]. The function g yields class invariants 
when D = 13 mod 24. In terms of our generic algorithm, we have n = 2, and 
for p = 2 mod 3 we get exactly two roots of ^g(A, j^) mod p, which differ only 
in sign. Thus H{X) = HD[g^]{X^) = i?D[ff](A)ffz)[ff](-A), and from this we 
easily obtain Hnlg^], and also Holg] if desired. 
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5 Computational Results 

This section provides performance data for the techniques developed above. We 
used AMD Phenom II 945 CPUs clocked at 3.0 GHz for our tests; the software 
was implemented using the gmp [22] and zn_poly [24] libraries, and compiled 
with gee [19]. 

To compute the class polynomial Holf], we require a bound on the size of 
its coefficients. Unfortunately, provably accurate bounds for functions / other 
than j are generally unavailable. As a heuristic, we take the bound B on the 
coefficients of Hd given by [30, Lcm. 8], divide log2 B by the asymptotic height 
factor c(/), and add a "safety margin" of 256 bits. Wc note that with the CM 
method, the correctness of the final result can be efficiently and unconditionally 
confirmed [5] , so we are generally happy to work with a heuristic bound. 

5.1 Class polynomial computations using the CRT method 

Our first set of tests measures the improvement relative to previous computa- 
tions with the CRT method. We used discriminants related to the construction 
of a large set of pairing-friendly elliptic curves, see [30, §8] for details. We re- 
constructed many of these curves, first using the Hilbert class polynomial Hd, 
and then using an alternative class polynomial Hd [/] . In each case we used the 
explicit CRT to compute Hd or HdI/] modulo a large prime q (170 to 256 bits). 

Table 1 gives results for four discriminants with \D\ » 10^", three of which 
appear in [30, Table 2] . Each column lists times for three class polynomial com- 
putations. First, we give the total time Ttot to compute Hd mod q, including 
the time Tenum spent enumerating Ell£)(Fp), for all the small primes p, using 
Algorithm 2 as it appears in §2.2. We then list the times T^num ^tot ob- 
tained when Algorithm 2 is modified to use gcd computations whenever it is 
advantageous to do so, as explained in §2.3. The gcd approach typically speeds 
up Algorithm 2 by a factor of 2 or more. 

For the third computation we selected a function / that yields class invariants 
for D, and computed HdI/] mod q. This polynomial can be used in place of Hd 
in the CM method (one extracts a root Xq of HdIJ] mod q, and then extracts 
a root of ^'/(a;o, J) mod q). For each function / we give a "size factor", which 
approximates the ratio of the total size of Hd to HdIJ] (over Z). In the first 
three examples this is just the height factor c(/), but in Example 4 it is 4c(/) 
because the prime 59 is ramified and we actually work with the square root of 
-ff£)[^59], as noted in §4.1, reducing both the height and degree by a factor of 2. 

We then list the speedup Ttot/^totl/l attributable to computing HdIJ] rather 
than Hd- Remarkably, in each case this speedup is about twice what one would 
expect from the height factor. This is explained by a particular feature of the 
CRT method: The cost of computing Hd mod p for small primes p varies signif- 
icantly, and, as explained in [30, §3], one can accelerate the CRT method with a 
careful choice of primes. When fewer small primes are needed, we choose those 
for which Step 1 of Algorithm 1 can be performed most quickly. 

The last line in Table 1 lists the total speedup T^ot/Tl^^lf] achieved. 
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Example 1 


Example 2 


Example 3 


Example 4 


\D\ 
h{D) 
riog2 B] 




13569850003 

20203 
2272564 

^^20203^ 


11039933587 

11280 
1359134 


12901800539 

54706 
5469776 

(327038 52) 


12042704347 

9788 
1207412 

(292447 3^2 432) 


Tenum (roots) 

Ttot 




6440 
19900 


10200 
23700 


10800 
52200 


21700 
42400 


T'enum (gcds) 
Ttot 




2510 
15900 


2140 
15500 


3440 
44700 


4780 
25300 


Function / 
Size factor 
TUf] 




36 
213 


^47 
24 
305 


A71 
36 
629 


120* 
191 


Speedup {TU/TU 
Speedup {Ttot/TU 


[/]) 75 
[/]) 93 


51 
78 


71 
83 


132 
222 



Table 1. Example class polynomial computations (times in CPU seconds) 



5.2 Comparison to the complex analytic method 

Our second set of tests compares the CRT approach to the complex analytic 
method. For each of the five discriminants hsted in Table 2 we computed class 
polynomials Holf] for the double Ty-quotient 1153,13 and the Weber f function, 
using both the CRT approach described here, and the implementation [14] of 
the complex analytic method as described in [12]. With the CRT we computed 
Hoif] both over Z and modulo a 256-bit prime q; for the complex analytic 
method these times are essentially the same. 







complex 


analytic 


CRT 




CRT mod q 




h{D) 


rD3,13 


f 


t03,13 


f 


It's, 13 


f 


6961631 


5000 


15 


5.4 


2.2 


1.0 


2.1 


1.0 


23512271 


10000 


106 


33 


10 


4.1 


9.8 


4.0 


98016239 


20000 


819 


262 


52 


22 


47 


22 


357116231 


40000 


6210 


1900 


248 


101 


213 


94 


2093236031 


100000 


91000 


27900 


2200 


870 


1800 


770 



Table 2. CRT vs. complex analytic (times in CPU seconds) 



Wc also tested a "worst case" scenario for the CRT approach: the discrim- 
inant D = —85702502803, for which the smallest non-inert prime is £1 = 109. 
Choosing the function most suitable to each method, the complex analytic 
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method computes ffu [1x1109,127] in 8310 seconds, while the CRT method com- 
putes Hd[Ai3i] in 7150 seconds. The CRT approach benefits from the attractive 
height factor of the Atkin functions, c{Aisi) = 33 versus 0(11)109,127) ~ 12.4, and 
the use of gcds in Algorithm 2. Without these improvements, the time to com- 
pute Hd with the CRT method is 1460000 seconds. The techniques presented 
here yield more than a 200-fold speedup in this example. 

5.3 A record-breaking CM construction 

To test the scalability of the CRT approach, we constructed an elliptic curve 
using \D\ = 1000000013079299 > 10^^, with hiD) = 10034174 > 10^. This 
yielded a curve — — 3x + c of prime order n over the prime field ¥q, where 

c = 12229445650235697471539531853482081746072487194452039355467804333684298579047; 
q = 28948022309329048855892746252171981646113288548904805961094058424256743169033; 
n = 28948022309329048855892746252171981646453570915825744424557433031688511408013. 

This curve was obtained by computing the square root of i?D[^7i] modulo q, a 
polynomial of degree h{D)/2 = 5017087. The height bound of 21533832 bits was 
achieved with 438709 small primes p, the largest of which was 53 bits in size. 
The class polynomial computation took slightly less than a week using 32 cores, 
approximately 200 days of CPU time. Extracting a root over Fg took 25 hours 
of CPU time using NTL [29]. 

We estimate that the size of ^/Hd\Ati\ is over 13 terabytes, and that the 
size of the Hilbert class polynomial Hu is nearly 2 petabytes. The size of 
^^HdIAyi] mod q, however, is under 200 megabytes, and less than 800 megabytes 
of memory (per core) were needed to compute it. 
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